Interplay of instabilities in mounded surface growth 
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We numerically study a one-dimensional conserved growth equation with competing linear 
(Ehrlich-Schwoebel) and nonlinear instabilities. As a control parameter is varied, this model ex- 
hibits a non-equilibrium phase transition between two mounded states, one of which exhibits slope 
selection and the other does not. The coarsening behavior of mounds in these two phases is studied 
in detail. In the absence of noise the steady-state configuration depends crucially on which of the 
two instabilities dominates the early time behavior. 
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The phenomenon of formation and coarsening of 
mounds in epitaxially grown thin films is a subject of 
much recent experimental 0, L| and theoretical [j. [j- UJ 
interest. Traditionally, the formation of mounds has been 
attributed to the presence of an Ehrlich-Schwoebel (ES) 
step-edge barrier 0,13 that hinders the downward motion 
of atoms across the edge of a step. This leads to an effec- 
tive "uphill" surface current 9] that has a destabilizing 
effect, leading to the formation of mounded structures. 
The ES mechanism is usually represented in continuum 
growth equations as a linear instability 9] that is con- 
trolled by higher-order nonlinear terms. In ES-type mod- 
els, the characteristics of the coarsening process, in which 
the typical lateral size of the mounds grows in time, de- 
pends on the details of the model. In these models, slope 
selection occurs (i.e. the slope of the mounds remains 
constant during coarsening) only if the "ES part" of the 
slope dependent surface current has one or more stable 
zeros as a function of the slope. 

In our earlier work [Tol Hl| on a class of spatially dis- 
crete, conserved, one-dimensional (Id) models of epitax- 
ial growth, we found a mechanism of mound formation 
and coarsening with slope selection that is different from 
the conventional ES mechanism. We studied the spa- 
tially discretized Lai-Das Sarma equation [12| of MBE 
growth and an atomistic model |13| that provides a dis- 
crete realization of the dynamics described by this equa- 
tion, and found the occurrence of a nonlinear instabil- 
ity |14| in which isolated pillars or grooves grow in time 
if their height or depth exceeds a "critical" value. When 
this instability is controlled by the introduction of an 
infinite number of higher order gradient nonlinearities, 
these models show, for a range of parameter values, the 
formation of mounds with a well-defined slope that re- 
mains constant during the coarsening process. These re- 
sults demonstrate that mound formation and power-law 
coarsening with slope selection can occur in the absence 
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of an ES instability. 

In most experimentally studied systems, however, it 
is believed that the ES step-edge barrier is present, al- 
though it may possibly be very weak. It is, therefore, 
important to understand how the behavior of our mod- 
els would be modified when the ES mechanism is in- 
corporated in their kinetics. To address this issue, we 
have studied, using numerical integration, a spatially dis- 
cretized Id growth equation in which a linear ES-type 
instability is present in conjunction with the nonlinear 
instability mentioned above. The main results of this 
study are as follows. Depending on the values of the 
parameters appearing in the model, we find two differ- 
ent mounded states. If the parameters are such that the 
nonlinear instability is the dominant one, then the behav- 
ior of the system is similar to that found in our earlier 
studies 0, ^3 : it exhibits the formation of triangular 
mounds and power-law coarsening with slope selection. 
If, on the other hand, the linear instability is dominant, 
then the system exhibits a different kind of mounded 
state in which the mounds have a cusp-like shape and 
they steepen during the coarsening process. We call 
these two mounded states "faceted" and "cusped", re- 
spectively. As the parameters are changed, the system 
undergoes a dynamical phase transition from one of these 
mounded states to the other. 

We study a spatially discretized Id version of the 
fourth order conserved growth equation, proposed in the 
context of MBE growth UQjl, in which the nonlinear in- 
stability 0| is controlled using a control function |ToL fll| 
of the form proposed by Politi and Villain |f| and a linear 
instability of the form proposed by Johnson et al. @ to 
represent the ES effect is also included. Thus, the equa- 
tion of motion of the interface height in appropriately 
non-dimcnsionalized form is written as: 

dhi/dt = -V 4 /i, + AV 2 (|V/i l | 2 /(l + ci|V/i,| 2 )) - 

V(V/i l /(l + c 2 |V/i l | 2 ))+ Vt (t), (1) 

where hi is the non-dimensionalized height variable at 
lattice site i, and V and V 2 are the lattice versions of 
the derivative and Laplacian operators, respectively, cal- 
culated using the nearest neighbors as outlined in our 
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earlier papers^], In Eq.QJ, ci, C2 and A are con- 
stants (model parameters), and rji(t) represents uncor- 
related random noise with zero mean and unit variance. 
Our results are based on numerical integration of this 
equation in Id, using a simple Eulcr scheme|14j in which 
the time evolution of the height variables is given by 

hi(t + At) = hi(t) + Ai{V 2 [-V 2 ^(i) + 
A(|V/i l (t)| 2 /(l + ci|V^(t)| 2 ))]- 
V(VAi(t)/(l + c 2 \Vhi(t)\ 2 ))} + VAl Vl (t). (2) 

In all our calculations, the value of the parameter C2 was 
held fixed at unity, and c\ and A were varied. Our results 
are obtained for different system sizes, 40 < L < 1000, 
with periodic boundary conditions. We do not find 
any significant L-dependence. We have used a time 
step At=0.01 for most of our studies and checked that 
very similar results are obtained with smaller values of 
At. The results obtained by using a more sophisti- 
cated algorithm 0] closely match with those of our 
Euler scheme for a small enough At. We have also 
checked that the basic features of our results remain un- 
altered if other boundary conditions (such as "fixed" and 
"zero-flux" boundary conditions considered in our earlier 
work 0,|n|) are used. 

We first describe the results obtained in the absence of 
the noise term. In this case, we find that if the parameter 
ci is sufficiently small, then the long time behavior of the 
interface depends on which of the instabilities dominates 
at early times. In order to characterize this, we consid- 
ered initial configurations with a single pillar of height ho 
on an otherwise flat interface, and studied the long-time 
behavior as a function of ho. We find that when ho is suf- 
ficiently small, so that the nonlinear instability is not ini- 
tiated, the linear instability dominates the time evolution 
and the resulting morphology is mounded without slope 
selection. If on the other hand ho is large enough to seed 
the nonlinear instability, then mounds with a "magic" 
slope result. This is illustrated in Fig. ^ which shows 
typical configurations obtained in the two cases: slope se- 
lected mounds ( "faceted" ) in situations where the height 
of the central pillar is greater than a critical value and 
mounds without slope selection ("cusp" like ) otherwise. 
Fig. also shows the dependence of the critical value of 
ho on A, the strength of the nonlincarity. We find that 
this critical value is approximately given by A/X, A ~ 21, 
with a weak dependence on the parameter c\. 

This "bistable" behavior is found only if the control 
parameter c\ is sufficiently small. Otherwise, the system 
evolves to the cusped morphology even if the initial con- 
figuration has a high central pillar. Thus, the faceted 
morphology is found in the zero-noise simulations only 
if ci < c 1 p (A) (the superscript sp denotes "spinodal", 
see below) and the initial configuration is sufficiently 
rough to seed the nonlinear instability. This behavior 
may be understood from a linear stability analysis. In 
both faceted and cusped growths, the finite-sized system 
evolves, at long times, to a profile with a single mound 
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FIG. 1: Dependence of the morphology on the initial con- 
ditions in the absence of noise. Main plot: the critical value 
of ho (see text) as a function of the nonlinearity parameter A 
for ci = 0.01. This critical value is approximately given by 
ho — 21/A (solid line). Upper inset: "faceted" morphology 
for h > ho. Lower inset: "cusped" morphology for h < ho. 



which is a fixed point of the noiseless dynamics. Exam- 
ples of such profiles (obtained from simulations with noise 
which slightly roughens the profiles) are shown in Figs. [21 
and |3 These fixed points can also be obtained by calcu- 
lating the hi for which <?;=0 for all i, where gi is the term 
multiplying At on the right-hand side of Eq.(J2Jl. The 
local stability of the faceted fixed point may be deter- 
mined from a calculation of the eigenvalues of the matrix 
Mij = dgi/dhj evaluated at the fixed point. We find 
that the largest eigenvalue of this matrix crosses zero 
at a "spinodal" value, c\ = c^ p (A) (see inset of Fig. [JJ, 
signaling an instability of the faceted profile. Thus, for 
< ci < c^(A), the dynamics of Eq.© without noise 
admits two locally stable invariant profiles: a cusped pro- 
file without slope selection, and a faceted one with slope 
selection. Depending on the initial state, the no-noise dy- 
namics takes the system to one of these two fixed points. 
For example, an initial state with one pillar on a flat back- 
ground is driven by the no-noise dynamics to the cusped 
fixed point if the height of the pillar is smaller than a 
critical value (mentioned earlier), and to the faceted one 
otherwise. The dependence of on the nonlinearity pa- 
rameter A is shown in Fig. EI Such a "spinodal line" does 
not exist for the cusped fixed point. 

We have studied in detail the process of coarsening 
of the mounds in the two different (faceted and cusped) 
growth modes. In these simulations, the initial configu- 
ration is obtained by setting hi = hoTi where n is a ran- 
dom number uniformly distributed between —1 and 1. If 
ci < Ci and ho is sufficiently large to initiate the nonlin- 
ear instability, the system evolves to a faceted structure; 
a cusped structure is obtained otherwise. The results 
reported here were obtained from averages over 20 such 
runs with different initial configurations. In the region 
of parameter space, < c < c^ p (A), where the faceted 
phase is locally stable, the mounds coarsen in time with 
the slope of the mounds remaining constant. In this case, 
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FIG. 2: The interface profiles at three different times (t= 
100, 4000 and 256000) in a run with noise starting from a flat 
state for an L — 400 sample with A = 4.0 and Ci = 0.05. A 
double log plot of the interface width W as a function of time 
t, averaged over 15 runs for L — 1000 samples in the presence 
of noise, and a power-law fit are shown in the inset. 




FIG. 4: Double log plots of the width W (full line) and 
the average of the maximum slope s m (dash-dotted line) as 
a function of time t, for faceted growth with slope selection 
in the absence of noise. This data are for A = 4.0, c\ — 0.01, 
and L = 800, averaged over 20 initial conditions. In this case 
the growth exponent j3 = 0.5 ± 0.01 (the dashed line is the 
best power-law fit) and the steepening exponent 9 — 0. 
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FIG. 3: Interface profiles at three different times (t= 200, 
4000 and 2560000) in a run starting from a flat state for an 
L = 400 sample with A = 4.0 and cj = 0.01. A double log 
plot of the interface width W as a function of time t, averaged 
over 15 runs for L = 1000 samples in the presence of noise, 
and power-law fits to the early and late-time data are shown 
in the inset. 




FIG. 5: Variation of width W (full line) and the average of 
the maximum slope, s m (dash-dotted line) as a function of 
time t on a log-log scale for cusped growth in the absence 
of noise (A = 4.0, ci = 0.05, and L = 800 averaged over 20 
initial conditions). The growth exponent is j3 = 0.54 ± 0.02, 
and the steepening exponent 8 — 0.2 ± 0.02. Best power-law 
fits are shown by the dashed lines. 



the average mound size R(t) is obviously proportional to 
the interface width W(t). As shown in FigQJ the inter- 
face width in this growth mode increases as a power-law 
with time in the long-time coarsening regime, while the 
slope of the mounds remains constant in time. The ex- 
ponent (3 that described this power-law coarsening be- 
havior, W(t) oc R{t) oc t 13 , is found to be close to 0.5. 
The coarsening of mounds in the cusped regime (i.e. for 
ci > c & i(X) and any initial configuration, and c\ < c\ p (A) 
and sufficiently smooth initial configurations) is qualita- 
tively different. In this growth mode, both the interface 
width and the average slope s(t) of the cusp-like mounds 
(as well as the maximum slope s m (t)) increases with time 
as a power law in the coarsening regime: W(t) oc t@ , 
and s(t) tx t e . This implies that the average mound size 
also increases with time as a power law: R(t) oc t n with 



n — (3 — 9. This behavior is illustrated in Fig. EJ The val- 
ues of the exponents are found to be (3 = 0.54 ± 0.02 and 
9 = 0.2 ± 0.02, implying that the coarsening exponent is 
n = 0.34 ±0.04. 

In the presence of noise, the steady state behavior is 
independent of the initial condition. When the control 
parameter c\ is sufficiently large, the nonlinear instabil- 
ity is completely suppressed and the route to mounding is 
similar to the linear instability dominated behavior men- 
tioned above, the configurations being slightly roughened 
versions of their noiseless counterpart. This is illustrated 
in Fig.|2]where the interface profiles in a typical run start- 
ing from a flat state for L = 400, A = 4.0 and ci = 0.05 
are shown at times t — 100 (early time regime), t = 4000 
(coarsening regime) and t = 256000 (steady state). The 
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FIG. 6: Phase diagram for a L= 200 system, where the crit- 
ical values of the control parameter ci are shown as functions 
of A. The critical value Cj P (ci r ) above which the faceted phase 
is unstable (metastable) is shown by squares (triangles). The 
dashed lines are guides to the eye. Left Inset: Zero crossing of 
the largest eigenvalue re of the stability matrix of the faceted 
fixed point as a function of ci (A = 4.0, L — 200). Right Inset: 
The probabilities Pi (circles) and P2 (inverted triangles) (see 
text) as functions of ci with A = 5.0 and L= 200. 



inset shows the interface width If as a function of t, 
averaged over 15 runs for L = 1000 samples. The aver- 
aged data show W oc t 13 , [3 = 0.5 ± 0.01 in the coarsen- 
ing regime. The average slope of the mounded interface 
grows as s(t) oc t e , with 9 = 0.18 ± 0.02 and hence the 
coarsening exponent is n = 0.32 ± 0.03. 

As the value of c\ is decreased below a critical value 
holding A fixed, the nonlinear instability dominates over 
the linear one. At early times in runs starting from a flat 
state, the interface is self affine and the interface width 
shows power law scaling with an exponent close to 0.37. 
As time progresses isolated pillars with height ho > h m i n , 
where h m i n is the minimum height of such a pillar above 
which the nonlinear instability is operative, make their 
appearance through random fluctuations. The time evo- 
lution of the interface beyond the point of occurrence 
of the instability is similar to that in the noiseless sit- 
uation. In Fig. 13 we present snapshots of an L = 400 
sample for A = 4.0 and c\ = 0.01 at different stages of 
growth: t = 200 (before the onset of the instability), 
t = 4000 (coarsening regime), and t — 2560000 (steady 
state). The inset shows a plot of the interface width as a 
function of time, obtained by averaging over 15 runs for 
L = 1000 samples. The averaged data show a power-law 
growth regime with exponent 0.37 ± 0.01 before the on- 
set of the instability and a second power-law coarsening 
regime with W oc , /3 = 0.49 ± 0.02. 

A dynamical phase transition at c\ = c^ r (A) < c 1 p (A) 
separates these two kinds of growth modes. To calcu- 
late cf(X), we start a system at the faceted fixed point 
and follow its evolution according to Eq.J2J for a long 
time (typically t — 10 4 ) to check whether it reaches a 
cusped steady state. By repeating this procedure many 
times (typically 100 runs), the probability, Pi(A,ci), of a 



transition to a cusped state is obtained. For fixed A, Pi 
increases rapidly from to 1 as c\ is increased above a 
critical value. Typical results for P% as a function of c\ 
for A = 5.0 are shown in the right inset of Fig. |SJ The 
value of c\ at which P\ = 0.5 provides an estimate of 
cf. Another estimate is obtained from a similar calcula- 
tion of P2(X,ci), the probability that a flat initial state 
evolves to a faceted steady state. As expected, Pi in- 
creases sharply from to 1 as c\ is decreased (see right 
inset of Fig. |SJ), and the value of C\ at which this prob- 
ability is 0.5 is slightly lower than the value at which 
P\ = 0.5. This difference reflects finite-time hysteresis 
effects. The value of cf is taken to be the average of 
these two estimates, and the difference between the two 
estimates provides a measure of the uncertainty in the 
determination of cf r . The phase boundary obtained this 
way is shown in Fig. |HJ along with the results for c^ p . 

The scaling behavior in the coarsening regime in the 
presence of noise is the same as that found in the noise- 
less case. The qualitative behavior in the faceted phase 
is similar to that found in our earlier work |lfjt lll| on 
this model without the ES term. The ES term, however, 
has an important effect: it changes the coarsening expo- 
nent from 0.33 to 0.5. A model very similar to the one 
considered here was studied by Torcini and Politi 
for parameter values deep in the cusped regime (small 
A, large c\). The mounded morphology we find in this 
regime is similar to that found in their study. Our results 
for the exponents 9 and n are slightly different from the 
values (9 = n = 1/4) reported by them. This is probably 
due to crossover effects - we have found values of n closer 
to 1/4 if smaller values of A and/or larger values of c\ are 
used. 

In summary we have shown that a nonlinear instability 
in our spatially discrete growth model in which the linear 
ES instability is also present may, for appropriate param- 
eter values, lead to mound formation with slope selection 
and power-law coarsening. This is qualitatively different 
from the behavior in the parameter regime where the ES 
instability dominates: the system exhibits mound forma- 
tion in this regime also, but there is no slope selection. 
The coarsening exponent has different values in these two 
regimes which are separated by a line of first-order dy- 
namical phase transitions. The ES part of the surface 
current in our model does not vanish for any non-zero 
value of the slope. Therefore, the slope selection we find 
in the regime where the nonlinear instability is dominant 
is qualitatively different from that in ES-type models and 
is a true example of nonlinear pattern formation. The 
noiseless version of our model exhibits an interesting de- 
pendence on initial conditions: the long-time behavior 
depends on whether the inhomogencitics in the initial 
configuration are sufficient to seed the nonlinear insta- 
bility. Both kinds of mounding behavior found in this 
study have been observed in experiments, and our model 
may be relevant in the development of an understanding 
of these experimental observations. However, more work 
is required to establish a connection between our results 
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and those of experiments. 
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